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We present a hydrodynamic model that captures the essence of granular 
dynamics in a vibrating bed. We carry out the linear stability analysis and 
uncover the instability mechanism that leads to the appearance of the con- 
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vective rolls via a supercritical bifurcation of a bouncing solution. We also 

o ' 

explicitly determine the onset of convection as a function of control parameters 
t>X)' and confirm our picture by numerical simulations of the continuum equations. 

9 ' PACS numbers: 47.20.-k, 46.10+z, 81.35. +k 

Granular materials in a container subjected to vertical vibrations display interesting nonlin- 
ear dynamical behaviors. Nothing really happens for T = Au 2 / g < 1, where A and uj 
are the frequency and the amplitude of the oscillations and g is the gravitational constant. 
But for 1 < T, the granular materials collectively move up and down, which we term the 
uniform bouncing @, until T reaches the critical value T c beyond which such a uniform 
bouncing motion becomes unstable and the permanent convective rolls develop inside the 
bulk. Recent studies have revealed further complexity of this problem for values of 

T much larger than one, where reverse convection || and bubble formation || have been 
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observed. Current efforts to understand the experiments of granular dynamics [0,|5|,§ have 
mostly focused on large scale Molecular Dynamics (MD) simulations. 0] While successful 
in reproducing convection cells and some of the experimental results, such studies have 
limitations in understanding the analytic structure of the instability mechanism and/or its 
subsequent dynamic evolution. There have been a handful of attempts in the past to derive 
continuum equations for granular dynamics notably by Jenkins and Savage for rapid granu- 
lar flow problems || and by Haff for vibrating beds ||, but these studies have been mostly 
confined either to simple cases of one dimensional oscillations in an infinite system where 
pressure inside the grains behave like a fluid |§, or to cases where explicit assumption has 
been made regarding the Gaussian nature in velocity distribution of grains ||, which have 



been shown to break down in a dense granular system ||10|| . There also has been a recent 
attempt by Bourzutchky and Miller JTI] who have utilized the Navier Stoke equation along 
the similar line of Haff and have reproduced numerically the convective rolls. However, 
we find it difficult to imagine that the hydrodynamic pressure term (pgz) exists to cancel 
the gravity term inside the granular materials that undergo vertical vibrations. 

The purpose of this Letter is two fold: we first propose a dynamic model that is simple 
enough to make progress in analytic studies, yet captures, in our opinion, the essence of 
granular dynamics of vibrating beds. Second, we demonstrate that the correct way of 
studying the convective instability is to carry out the stability analysis around the bouncing 
solution and explicitly determine the onset of convection as a function of external parameters. 
We will also present numerical results to confirm our predictions. 

Equations of motion: Our starting point is the recognition that the most fundamental 
aspect of the vibrating bed, apart from the obvious fixed bed solution with no external 
driving, is the existence of a uniform bouncing of a collection of particles. Such a bouncing 
solution can be either a solid block inside the bed or a fluidized state with a slightly ex- 
panded volume yet with no internal degrees of freedom. This assumption is consistent with 
observations in MD M where surface fluidization rapidly spreads out into bulk regions when 



2 



surface fluidization is suppressed. In such a case, the bouncing solution can be represented 
by a motion of a ball on a vibrating platform. For small T, no exotic motion such as chaotic 
motion is expected to occur for such a ball [Hfl. We further assume that the restitution 



constant of the collection of particles is zero to represent the relaxation of inside collection 
of particles. In such a case, the relative position of the ball with respect to the bottom plate, 
A(t), is given by: 

A(t) = r(sint - sin*) + T cos t (t — t ) - -(t - t ) 2 (1) 

in the unit of g — ou — 1, where the ball starts to bounce at to on the bottom plate, whose 
position at time t is given by T sin t in the experimental frame. The bouncing solution is 
then described by the relative speed between the plate and the ball: V re i = dA(t)/dt. Since 
A(t) cannot be negative, the ball launched upward on the plate at to falls back to the plate at 
ti(i.e,A(ti) = 0) and stays there until t = to + 2ir from our assumption of the zero restitution 
constant. The ball is then relaunched and obeys (1) again. For later use, we determine 
{t = sin _1 (l/r),ti) for different values of T. For example, (t„,ti) = (1.181,2.88225) for 
r = 1.1 and (t , ti) = (0.524,5.18) for T = 2.0. When we expand A(t) around t we 
obtain A(t) ~ (r/6) cos(t )(t — to) 3 > 0, where cos(t ) > from the launching condition 
d 2 V re i/dt 2 > 0. Hence, there is no solution of A(t) = around t = to except for t = t . 
Therefore, the bouncing motion starts from the finite t x — t . One can now readily derive 
the equation of motion for the vertical coordinate z for the bouncing motion of a granular 
block: z = (—1 + r sin t)6 l (— 1 + Tsint) where the 8(x) = 1 for x > and 8(x) = for 
otherwise. 

In order to describe the motion of a granular block in the presence of internal degrees 
of freedom such as rotation and/or translation, we define two coarse-grained dynamical 
variables: the density p(r, t) and the velocity v(r,t) of the granular system. In a box fixed 
frame, p and v then should satisfy the continuity and the momentum equation: 

d t p + V • (pv) = (2) 
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d t v + (v • V)v = z(T sin t - 1 - A) - - VP + - [V 2 v + X V( V • v)] (3) 

p it 

where 5 is the unit vector in the vertical direction and A is a Lagrange multiplier. A = for 
free motion and A = T sin t — 1 for stationary state. Note that the first term in (|3]) is due 
to the uniform bouncing and the third term is the energy dissipation effectively represented 
by the Reynolds number R and the bulk viscosity x- Now, the exact form of the pressure P 
in (|3|) is unknown for granular materials. Unlike fluid, for granular materials in a container 
supported by the side walls, the pressure inside the bulk seems to saturate PUTS]. In such 
a case, the only contribution to the granular pressure would result from the hard sphere 
repulsion which might be effectively represented by the Van der Waals equation: 

where T k-< v 2 > is the granular temperature and b is a constant of order unity. 
Note that eqs.(fj) and (H) are precisely the compressible Navier-Stokes equation with two 
modifications: first the hydrodynamic pressure term is absent, which is replaced by the Van 
der Waals form (H), and second, the gravity term thus survives in the vibrating bed and has 
been effectively modified to g — Auj 2 sin t in the physical unit. Notice that the term Tsint 
appears since we have used the box- fixed frame. We now analyze eqs.fl^) and (|3|). 

Linear stability analysis: (a) Fixed bed solution: Fixed bed is a container filled with 
grains with no external driving. In this case, the contact force balances out the gravity and 
the net force acting on each grain is zero. So, we use A = F sin t — 1 in (|3|). In this case, the 
solution with constant density p = po and zero speed v = is stable. 

(b) Linear stability of a uniform bouncing solution: In order to discuss the stability of the 
uniform bouncing solution, p = p and v = (0, 0, V re i(t)), against fluctuations, we set p = p + 
pi, and decompose the velocity into the vertical and horizontal components, = (vj^, Wl) 
with = vj_ t i and w = V re i(t) +wl- We then substitute these into dynamic equations (0) 
and (^|) and introduce a new coordinate to simplify the problem, £ = z — J 1 V re i(t')dt'. Upon 
linearization, we obtain the following equation for the perturbed density pi- 



[d 2 - ^V 2 «9 f - T e V\}p L = (5) 

where R = R/(l + x)- We now solve ([5]) in 2 dimension under the no current boundary 
condition at the plate and at z — oo, namely: Pl = at x — 0, L, and z = 0, oo 
where L is the dimensionless size of the box. To satisfy these boundary conditions, we set: 

p L (x, y, z, t) = PL, q ,m(t) sin[7rmx] sin[g(£ - S(t))] (6) 

where 7r = ir/L, m is an integer, and S(t) = —A(t). We notice that the spectrum is discrete 
for x direction but continuous along z direction. We now substitute @ into fl5|) and utilize 
the fact, t = r + to with to — sin -1 1/r. After some algebra, we obtain the following second 
order ordinary differential equation for the amplitude p q (t) = pL,q,m(t): 

p q + B(q)p q + iC(q)p q + D(q)p q + iE{q)p q = %L q {r)p q + M q {r)p q + iN q (r)p q (7) 

where 

B(q) = R-\7i 2 m 2 + q 2 ), C(q) = 2qVf 2 ~^T (8) 

D(q) = T e (q 2 + 7t 2 m 2 ) - ^q 2 T 2 + q 2 , E(q) = -q + ^/Y^lR^q^m 2 + q 2 ) (9) 

and the time dependent inhomogeneous terms are: L q (r) = 2q[r + \/T 2 — 1 cost — sinr], 
M q (r) = -2q 2 (T 2 - l)cosr + 2q 2 ^/T 2 - 1 sinr + ^^cos(2r) - q 2 ^T 2 - lsin(2r) - 
2q 2 Vr 2 — lr (1 — cos r) — 2q 2 r sin r + q 2 r 2 and N q (r) = —q\/T 2 — lsin r — q cos r + Rr 1 q(q 2 + 
7r 2 m 2 )(r + vT 2 — 1 cost — sinr). 

Note that the equation (^) is valid only between r = 2mr and r = r + 2n7r with 
ro = t\ — to, during which grains are launched from the plate by vibrations and then undergo 
free fall. Except for this region, it is easy to show S(t) = and C(q) = E(q) = L q (r) = 
M q (r) = N q (r) = with D(q) -> D (q) = T e (q 2 + vr 2 m 2 ). 

The rest of the paper is devoted to discuss the condition for the linear stability of (0) 
and numerically test the validity of such approximations. We may be able to obtain an 
explicit solution of Eq. (0) with the aid of assumption that the most unstable mode is only 
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a relevant mode. The condition for instability from this treatment, however, is complicated 
and time dependent. In addition, this condition is not adequate for our purpose, since we 
are interested in the behavior in time longer than one vibrating oscillation. Therefore, we 
may replace L q (r) (and for M q (r) and N q (r) as well) by its average value over a flying time 
< L q (r) >, namely L q (r) ~< L q >= i J T ° drLqij) and so on for M q {j) and N q (r). Eq.(^) 
is then reduced to a second order ordinary differential equation with constant coefficients. 
Assuming p q ~ e at , it becomes easy to obtain the eigenvalues a for the flying motion as 

B + iC , 1 



a ± = — ± - VOB + iCf - A(D + iE) (10) 

where C = C(q)- < L q >, D = D(q)- < M q > and E = E(q)- < N q >. The relevant 
branch is only cr + . On the other hand, the eigenvalues are reduced to a± = —B/2 ± 



\/B 2 — ADq/2 for stationary states. 

The averaged instability condition over one oscillation cycle is then the average of Re[a] > 
0. For this purpose, we introduce a function: 

**ff(Q) = M(E - ^) 2 - B\D + ^)} + (2tt - r )(-B 2 D ). (11) 

where the first term is the instability condition for ( PH| ) multiplied by the time period, To, 
in which particles can move freely ||14|| , while the second term is that with S(t) = 0. If the 
function cr e ff(q) > for any q, it signals the instability of the uniform bouncing solution. 
For finite system, cr e //(0) = —27i 7 T e /(R 2 L 6 ) < 0. Thus, the convection will disappear for 
infinite systems, which agrees with MD simulations Jl]|||15|]. Equivalently, convection also 
disappears in the limit of large R, i.e. either the particles are too smooth or the kinetic energy 
is too small to provide the necessary driving force among grains. The set of parameters that 
corresponds to physical situations might be :R ~ 2, T e ~ 3 and L = 10, because (i) the linear 
size of the box L r = Lg/uj 2 ~ 0.6[cm] for u ~ 20 [.Hz], (ii) T ~ / TO V? el (r)dT ~ 3, (iii) the 
kinetic viscosity for granular fluid is evaluated by u s ~ 5 x 10~ 3 [m 2 /s] Q and the definition 



of R = UL r /v s ~ 2 with the aid of the characteristic velocity U ~ yV 2 el g/u ~ 10cm/ sec 
in the physical unit. But for pure numerical reasons, we choose R = T e = 10 and L = 10. 



For this set of parameters, we first solve A(t) = numerically to determine t%, and then 
compute a e ff(q) as a function of q. As demonstrated in Fig.l, a e ff(q) is convex and thus has 
a maximum, a m , at a particular value of q. For r ~ 1, a m < and thus cr e //(g) < for all 
g and the bouncing solution is stable. But as we increase T further to the critical value, T c , 
a m moves upward crossing zero and becomes positive, in which case cr e ff{q) > for a band 
of q. In this case, the bouncing solution becomes unstable and we expect the convective 
rolls to appear. The onset of convection is then determined by setting, cr m (r c ) = 0. For 
L = 10, R = T e = 10, we find T c ~ 1.12 and the selcted wave number is about q c = 0.22. 
The most unstable wave number q m gradually shifts with the increase of V. We now check 
the validity of our picture by numerical simulations. 

Numerical Results: We have solved (2)-(4) numerically in two dimension with no slip 
boundary conditions at the side walls as well as at the top and the bottom plates. Note 
that the top plate suppress complicated surface motion of vibrating beds and allow us to 
use the simplified picture. Since the granular fluid is confined in a box, we do not introduce 
A explicitly in the simulations. As a result, S(t) ~ after a grain lands on a plate in the 
average bouncing state. The absence of A and the presence of the top wall is expected to 
cause the appearance of the bouncing solution for T < 1 in contrast with the real situation. 
But since the linearized eq.(^) with S(t) =0 is identical to that with non-zero A, omitting A 
would not change the essence of the dynamics. In the same spirit, we have ignored \ and b 
in our simulations. Our simulation results are presented in Fig. 2 for two different values of 
T, T < r c and T > r c . In the former case, the bouncing solution is expected to appear inside 
the bed and the density and the velocity at a given point oscillates with the same frequency 
of the vibration. (Fig. 2a) Upon increasing F further to T = 1.2, which is beyond the predicted 
r c = 1.12 determined by (11), we find that the bouncing solution has disappeared and the 
permanent convective rolls have developed inside the bulk (Fig.2b). The wavelength of the 
most unstable mode by the linear stability analysis is about q m ~ 0.4, which is not far from 
the actual wavelength of the convective rolls: q = 2n/\ = 2ix/L ph 0.6. 
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In passing, we briefly mentioned the difference between the granular beds and the water 



beds. The latter is shown to exhibit the Faraday instability on the air-water interface, flc 
The crucial difference between these two systems lie in the pressure term: for the water bed, 
since the water is an incompressible fluid, the hydrodynamic pressure term pgz precisely 
cancels the gravity term in the fluid equation, thus suppressing the motion inside the bulk, 
while the absence of the hydrodynamic pressure term produces the convective instability 
in the bulk for the granular beds. We will present the details of our analysis including the 
weakly nonlinear analysis elsewhere, which will highlight the differences between the two. 
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Figure Captions 



Fig.l. The effective growth rate cr e ff(q) as a function of the wave number q for T = 
1.05 (diamond), for which a e ff(q) < for all values of q, while for r = 1.2 > r c = 1.12, 
a e ff(q) becomes positive for a band of q(square). T c is determined by the condition that 
the maximum of <J e ff(q) becomes zero at T c . (cross) The parameters used are: T e = R = 10 
and L = 10. 

Fig. 2 (a) A bouncing solution. The speed v z at a given point is plotted as a function time 
for r = 0.9. (b) For T = 1.2 > T c = 1.12, the bouncing solution becomes unstable and 
the permanent convective rolls appear inside the box. The arrows are the velocity vectors 
pointing upward. The parameters used in simulations for (a) and (b) are the same as those 
in Fig.l. 
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